function [optX,variance ] = meanVar( prices, dt, Rp, v)
%CVAROPTIMIZE 
%   Function find optimal portfolio with mean variance constraints
%   using quadprog
%                      Inputs:
% prices: Historical prices with J senarios and n (rows) stocks
% Note: J here is having dt more senarios, so real J = J - dt;
% Note: No risk free asset will be included since this is for mean-var
%       efficient frontier
%
% dt: optimization period   (default 10 days)
% Rp: required expected portfolio return
% v:  value constraint

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% input test
if nargin < 2
    dt = 10; Rp = 0.01; v = 0.2;
elseif nargin < 3
    Rp = 0.01; v = 0.2;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function Starts

[n, J] = size(prices);
J = J - dt; % get back the index to be correct, since prices includes
            % dt more days of data for the optimization period

%rf_v = rf*ones(1,J);

Rij = prices(:,1+dt:J+dt)./prices(:,1:J) - 1;
%Rij = [rf_v; Rij];
Er = mean(Rij, 2);
sigma = cov(Rij');

[n, J] = size(Rij); % update the size of n

% Now set up for using 'quadprog'

f = [];

% Set up for Aeq
% constraint (29)
Beq1 = 1;
Aeq1 = ones(1, n);

% constraint (30)
Beq2 = Rp;
Aeq2 = Er';

Aeq = [Aeq1; Aeq2];
Beq = [Beq1; Beq2];

LB = zeros(n,1);
UB = v*ones(n,1);

% Before trying to optimize, check if expected return Rp is achievable.
max_r = max(Er);
assert(Rp <= max_r, sprintf('Max expected return is %f, cannot reach\n'...
    ,max_r));

% Solve the optimal problem by 'quadprog'
[optX,fval,exitflag] = quadprog(sigma,f,[],[],Aeq, Beq, LB, UB);
assert(exitflag == 1, sprintf('No feasible solution with flag: %d', exitflag));
variance = 2*fval;


end

